{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the analysis from 10/27, let's fix the errors that might be causing the duplicate variables. Also note that some of these results might be a bit different because today I found an error in the pedigree for trio 10215_trio3, so that's fixed now. I re-ran all SNV tools for that particular trio, and redid the ensembles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Within tools stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GATK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "suffix=hiConfDeNovo.vcf \n",
    "cd /data/NCR_SBRB/simplex/gatk_refine\n",
    "rm interesting_snvs_${suffix}.txt\n",
    "# figure out all family IDs\n",
    "ls -1 *_trio1_${suffix} > famids.txt\n",
    "sed -i -e \"s/_trio1_${suffix}//g\" famids.txt\n",
    "\n",
    "# for each family ID\n",
    "while read fam; do\n",
    "  # figure out how many trios we have\n",
    "   ntrios=`ls -1 ${fam}_trio?_${suffix} | wc -l`;\n",
    "   ntrios=$(($ntrios))\n",
    "   # if we have more than one (assuming the first one is affected)\n",
    "   if [ $ntrios -gt 1 ]; then\n",
    "      # get all SNVs in the affected trio in the family\n",
    "      cut -f 1,2 ${fam}_trio1_${suffix} | grep -v '#' - | sort | uniq > ${fam}_possible_snvs_${suffix}.txt;\n",
    "      # combine the vcf files of all unnafected trios\n",
    "      cat ${fam}_trio[2..$ntrios]_${suffix} > ${fam}_control_snvs_${suffix}.txt;\n",
    "      # for each possible SNV in affected trio, mark it as interesting if it's not\n",
    "      # in the unnafected trios\n",
    "      while read snv; do\n",
    "         if ! grep -q \"$snv\" ${fam}_control_snvs_${suffix}.txt; then\n",
    "            echo $snv >> interesting_snvs_${suffix}.txt;\n",
    "         fi;\n",
    "      done < ${fam}_possible_snvs_${suffix}.txt;\n",
    "   fi;\n",
    "done < famids.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "def do_perms(snps, noften, nboot=10000, npicks=20):\n",
    "    success = 0\n",
    "    for i in range(nboot):\n",
    "        picks = np.random.choice(snps, npicks, replace = True)\n",
    "        counts = stats.itemfreq(picks)\n",
    "        nmax = np.max(counts[:,1])\n",
    "        if (nmax >= noften):\n",
    "            success += 1\n",
    "    return(success/float(nboot))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum frequency: 2 in 3388 unique snps\n",
      "[['chr1 40229504' 2]\n",
      " ['chr17 72762902' 2]\n",
      " ['chr18 15271172' 2]\n",
      " ['chr2 89105006' 2]\n",
      " ['chr20 29637674' 2]\n",
      " ['chr20 29637691' 2]\n",
      " ['chr6 29968761' 2]\n",
      " ['chr6 32548712' 2]\n",
      " ['chr6 32557647' 2]\n",
      " ['chr7 154002420' 2]]\n"
     ]
    }
   ],
   "source": [
    "snps_gatk = pd.read_table('/data/NCR_SBRB/simplex/gatk_refine/interesting_snvs_hiConfDeNovo.vcf.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_gatk = stats.itemfreq(snps_gatk['snp'])\n",
    "my_max = np.max(counts_gatk[:, 1])\n",
    "print 'Maximum frequency: %d in %d unique snps' % (my_max, len(counts_gatk))\n",
    "print counts_gatk[counts_gatk[:, 1] == my_max]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0524"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "do_perms(counts_gatk[:, 0], 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The issue we had was not in GATK, but it's good that the changes we made did't affect the results. They're still not great, but maybe close?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Triodenovo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "suffix=denovo_v2.vcf \n",
    "cd /data/NCR_SBRB/simplex/triodenovo\n",
    "rm interesting_snvs_${suffix}.txt\n",
    "# figure out all family IDs\n",
    "ls -1 *_trio1_${suffix} > famids.txt\n",
    "sed -i -e \"s/_trio1_${suffix}//g\" famids.txt\n",
    "\n",
    "# for each family ID\n",
    "while read fam; do\n",
    "  # figure out how many trios we have\n",
    "   ntrios=`ls -1 ${fam}_trio?_${suffix} | wc -l`;\n",
    "   ntrios=$(($ntrios))\n",
    "   # if we have more than one (assuming the first one is affected)\n",
    "   if [ $ntrios -gt 1 ]; then\n",
    "      # get all SNVs in the affected trio in the family\n",
    "      cut -f 1,2 ${fam}_trio1_${suffix} | grep -v '#' - | sort | uniq > ${fam}_possible_snvs_${suffix}.txt;\n",
    "      # combine the vcf files of all unnafected trios\n",
    "      cat ${fam}_trio[2..$ntrios]_${suffix} > ${fam}_control_snvs_${suffix}.txt;\n",
    "      # for each possible SNV in affected trio, mark it as interesting if it's not\n",
    "      # in the unnafected trios\n",
    "      while read snv; do\n",
    "         if ! grep -q \"$snv\" ${fam}_control_snvs_${suffix}.txt; then\n",
    "            echo $snv >> interesting_snvs_${suffix}.txt;\n",
    "         fi;\n",
    "      done < ${fam}_possible_snvs_${suffix}.txt;\n",
    "   fi;\n",
    "done < famids.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum frequency: 5 in 18802 unique snps\n",
      "[['chr1 16914580' 5]\n",
      " ['chr10 49319323' 5]\n",
      " ['chr12 69667893' 5]\n",
      " ['chrX 102973509' 5]\n",
      " ['chrX 118751076' 5]\n",
      " ['chrX 13397236' 5]\n",
      " ['chrX 135307049' 5]\n",
      " ['chrX 14599572' 5]\n",
      " ['chrX 152610294' 5]\n",
      " ['chrX 153008911' 5]\n",
      " ['chrX 153880181' 5]\n",
      " ['chrX 153904473' 5]\n",
      " ['chrX 154774663' 5]\n",
      " ['chrX 16876980' 5]\n",
      " ['chrX 38262808' 5]\n",
      " ['chrX 43601142' 5]\n",
      " ['chrX 96502650' 5]]\n"
     ]
    }
   ],
   "source": [
    "snps_tdn = pd.read_table('/data/NCR_SBRB/simplex/triodenovo/interesting_snvs_denovo_v2.vcf.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_tdn = stats.itemfreq(snps_tdn['snp'])\n",
    "my_max = np.max(counts_tdn[:, 1])\n",
    "print 'Maximum frequency: %d in %d unique snps' % (my_max, len(counts_tdn))\n",
    "print counts_tdn[counts_tdn[:, 1] == my_max]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results here didn't change either... Do we really have 5 trios with each of these SNPs?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/data/NCR_SBRB/simplex/triodenovo/10094_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/10128_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/10153_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/10448_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/855_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/10153_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/10406_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/1892_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/1893_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/1976_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/10033_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/10153_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/10164_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/10178_trio1_denovo_v2.vcf\n",
      "/data/NCR_SBRB/simplex/triodenovo/1976_trio1_denovo_v2.vcf\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "grep -l \"16914580\" /data/NCR_SBRB/simplex/triodenovo/*trio1_denovo_v2.vcf\n",
    "grep -l \"49319323\" /data/NCR_SBRB/simplex/triodenovo/*trio1_denovo_v2.vcf\n",
    "grep -l \"69667893\" /data/NCR_SBRB/simplex/triodenovo/*trio1_denovo_v2.vcf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apparently yes. Not always the same families, but always 5 different ones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DenovoGear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "suffix=dnm.vcf \n",
    "cd /data/NCR_SBRB/simplex/dng\n",
    "rm interesting_snvs_${suffix}.txt\n",
    "# figure out all family IDs\n",
    "ls -1 *_trio1_${suffix} > famids.txt\n",
    "sed -i -e \"s/_trio1_${suffix}//g\" famids.txt\n",
    "\n",
    "# for each family ID\n",
    "while read fam; do\n",
    "  # figure out how many trios we have\n",
    "   ntrios=`ls -1 ${fam}_trio?_${suffix} | wc -l`;\n",
    "   ntrios=$(($ntrios))\n",
    "   # if we have more than one (assuming the first one is affected)\n",
    "   if [ $ntrios -gt 1 ]; then\n",
    "      # get all SNVs in the affected trio in the family\n",
    "      cut -f 1,2 ${fam}_trio1_${suffix} | grep -v '#' - | sort | uniq > ${fam}_possible_snvs_${suffix}.txt;\n",
    "      # combine the vcf files of all unnafected trios\n",
    "      cat ${fam}_trio[2..$ntrios]_${suffix} > ${fam}_control_snvs_${suffix}.txt;\n",
    "      # for each possible SNV in affected trio, mark it as interesting if it's not\n",
    "      # in the unnafected trios\n",
    "      while read snv; do\n",
    "         if ! grep -q \"$snv\" ${fam}_control_snvs_${suffix}.txt; then\n",
    "            echo $snv >> interesting_snvs_${suffix}.txt;\n",
    "         fi;\n",
    "      done < ${fam}_possible_snvs_${suffix}.txt;\n",
    "   fi;\n",
    "done < famids.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum frequency: 6 in 20901 unique snps\n",
      "[['chr2 89072915' 6]\n",
      " ['chr8 97156476' 6]]\n"
     ]
    }
   ],
   "source": [
    "snps_dng = pd.read_table('/data/NCR_SBRB/simplex/dng/interesting_snvs_dnm.vcf.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_dng = stats.itemfreq(snps_dng['snp'])\n",
    "my_max = np.max(counts_dng[:, 1])\n",
    "print 'Maximum frequency: %d in %d unique snps' % (my_max, len(counts_dng))\n",
    "print counts_dng[counts_dng[:, 1] == my_max]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, 6 families is better than 5 (or 2). But the other tools didn't pick that one up. Still worth analyzing it..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/data/NCR_SBRB/simplex/dng/10128_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/10131_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/10164_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/10173_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/10182_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/1893_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/10164_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/10178_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/10197_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/10406_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/1892_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/1976_trio1_dnm.vcf\n",
      "/data/NCR_SBRB/simplex/dng/855_trio1_dnm.vcf\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "grep -l \"89072915\" /data/NCR_SBRB/simplex/dng/*trio1_dnm.vcf\n",
    "grep -l \"97156476\" /data/NCR_SBRB/simplex/dng/*trio1_dnm.vcf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, all 6 different families. There's also a 7th family, but likely there the variable didn't get selected because it was also in the unaffected sib."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Within group stats\n",
    "\n",
    "The approach here is to calculate the best stats in all ADHD samples, and see what's the best we can do for a specific variable in non-affected siblings. For comparison, we can do it the other way around as well.\n",
    "\n",
    "## Triodenovo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "cd /data/NCR_SBRB/simplex/triodenovo/\n",
    "\n",
    "# concatenate all affected trios and extract the SNPs\n",
    "rm affected_snvs.txt unaffected_snvs.txt;\n",
    "for f in `ls *_trio1_denovo_v2.vcf`; do\n",
    "    cat $f | grep -v '#' - | awk 'BEGIN {FS=\"\\t\"; OFS=\":\"}; {print $1, $2}' - | sort | uniq >> affected_snvs.txt;\n",
    "done\n",
    "# need to be careful not to double count unaffected SNPs\n",
    "while read fam; do\n",
    "    if [ -e ${fam}_trio2_denovo_v2.vcf ]; then\n",
    "        cat ${fam}_trio[2..4]_denovo_v2.vcf | grep -v '#' - | awk 'BEGIN {FS=\"\\t\"; OFS=\":\"}; {print $1, $2}' - | sort | uniq >> unaffected_snvs.txt;\n",
    "    fi;\n",
    "done < famids.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum frequency: 14 in 21 affected families\n",
      "[['chrX:41093413' 14]]\n",
      "Maximum frequency: 11 in 29 unaffected families\n",
      "[['chrX:100098355' 11]\n",
      " ['chrX:150573743' 11]\n",
      " ['chrX:150575385' 11]\n",
      " ['chrX:27479339' 11]\n",
      " ['chrX:9686187' 11]]\n"
     ]
    }
   ],
   "source": [
    "naff = 21\n",
    "nunaff = 29\n",
    "aff_tdn = pd.read_table('/data/NCR_SBRB/simplex/triodenovo/affected_snvs.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_aff_tdn = stats.itemfreq(aff_tdn['snp'])\n",
    "my_max = np.max(counts_aff_tdn[:, 1])\n",
    "print 'Maximum frequency: %d in %d affected families' % (my_max, naff)\n",
    "print counts_aff_tdn[counts_aff_tdn[:, 1] == my_max]\n",
    "\n",
    "unaff_tdn = pd.read_table('/data/NCR_SBRB/simplex/triodenovo/unaffected_snvs.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_unaff_tdn = stats.itemfreq(unaff_tdn['snp'])\n",
    "my_max = np.max(counts_unaff_tdn[:, 1])\n",
    "print 'Maximum frequency: %d in %d unaffected families' % (my_max, nunaff)\n",
    "print counts_unaff_tdn[counts_unaff_tdn[:, 1] == my_max]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, we found a mutation in 14 out of the 21 affected families, but we can also find mutations happening in 11 out of 29 unaffected families. So, 14 is not terribly impressive. But let's see how often this one mutation happens in unaffected famlies?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([9], dtype=object)"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "counts_unaff_tdn[counts_unaff_tdn[:, 0] == 'chrX:41093413', 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can potentially assign some significance to that... but what if we strict ourselves to autossomal chromosomes?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum frequency: 9 in 21 affected families\n",
      "[['chr13:42894052' 9]\n",
      " ['chr3:39556476' 9]\n",
      " ['chrM:3012' 9]]\n",
      "Maximum frequency: 8 in 29 unaffected families\n",
      "[['chr1:7889972' 8]\n",
      " ['chrM:3012' 8]]\n",
      "Counts chr13:42894052 shows up in unaffected: 4\n"
     ]
    }
   ],
   "source": [
    "aff_tdn = pd.read_table('/data/NCR_SBRB/simplex/triodenovo/affected_snvs.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_aff_tdn = stats.itemfreq(aff_tdn['snp'])\n",
    "unaff_tdn = pd.read_table('/data/NCR_SBRB/simplex/triodenovo/unaffected_snvs.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_unaff_tdn = stats.itemfreq(unaff_tdn['snp'])\n",
    "\n",
    "keep_me = [s for s, snp in enumerate(counts_aff_tdn[:, 0]) if snp.find('chrX') < 0]\n",
    "counts_aff_tdn = counts_aff_tdn[keep_me, ]\n",
    "keep_me = [s for s, snp in enumerate(counts_unaff_tdn[:, 0]) if snp.find('chrX') < 0]\n",
    "counts_unaff_tdn = counts_unaff_tdn[keep_me, ]\n",
    "\n",
    "my_max = np.max(counts_aff_tdn[:, 1])\n",
    "max_str = counts_aff_tdn[counts_aff_tdn[:, 1] == my_max, 0][0]\n",
    "print 'Maximum frequency: %d in %d affected families' % (my_max, naff)\n",
    "print counts_aff_tdn[counts_aff_tdn[:, 1] == my_max]\n",
    "my_max = np.max(counts_unaff_tdn[:, 1])\n",
    "print 'Maximum frequency: %d in %d unaffected families' % (my_max, nunaff)\n",
    "print counts_unaff_tdn[counts_unaff_tdn[:, 1] == my_max]\n",
    "print 'Counts %s shows up in unaffected: %d' % (max_str, counts_unaff_tdn[counts_unaff_tdn[:, 0] == max_str, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DenovoGear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "cd /data/NCR_SBRB/simplex/dng/\n",
    "\n",
    "# concatenate all affected trios and extract the SNPs\n",
    "rm affected_snvs.txt unaffected_snvs.txt;\n",
    "for f in `ls *_trio1_dnm.vcf`; do\n",
    "    cat $f | grep -v '#' - | awk 'BEGIN {FS=\"\\t\"; OFS=\":\"}; {print $1, $2}' - | sort | uniq >> affected_snvs.txt;\n",
    "done\n",
    "for f in `ls *_trio[2..4]_dnm.vcf`; do\n",
    "    cat $f | grep -v '#' - | awk 'BEGIN {FS=\"\\t\"; OFS=\":\"}; {print $1, $2}' - | sort | uniq >> unaffected_snvs.txt;\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum frequency: 8 in 21 affected families\n",
      "[['chr17:21319860' 8]\n",
      " ['chrX:8432783' 8]]\n",
      "Maximum frequency: 7 in 29 unaffected families\n",
      "[['chr4:113190251' 7]\n",
      " ['chrX:13779124' 7]\n",
      " ['chrX:152772473' 7]\n",
      " ['chrX:8432783' 7]\n",
      " ['chrY:9967496' 7]]\n"
     ]
    }
   ],
   "source": [
    "naff = 21\n",
    "nunaff = 29\n",
    "aff_dng = pd.read_table('/data/NCR_SBRB/simplex/dng/affected_snvs.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_aff_dng = stats.itemfreq(aff_dng['snp'])\n",
    "my_max = np.max(counts_aff_dng[:, 1])\n",
    "print 'Maximum frequency: %d in %d affected families' % (my_max, naff)\n",
    "print counts_aff_dng[counts_aff_dng[:, 1] == my_max]\n",
    "\n",
    "unaff_dng = pd.read_table('/data/NCR_SBRB/simplex/dng/unaffected_snvs.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_unaff_dng = stats.itemfreq(unaff_dng['snp'])\n",
    "my_max = np.max(counts_unaff_dng[:, 1])\n",
    "print 'Maximum frequency: %d in %d unaffected families' % (my_max, nunaff)\n",
    "print counts_unaff_dng[counts_unaff_dng[:, 1] == my_max]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum frequency: 8 in 21 affected families\n",
      "[['chr17:21319860' 8]]\n",
      "Maximum frequency: 7 in 29 unaffected families\n",
      "[['chr4:113190251' 7]\n",
      " ['chrY:9967496' 7]]\n",
      "Counts chr17:21319860 shows up in unaffected: 6\n"
     ]
    }
   ],
   "source": [
    "aff_dng = pd.read_table('/data/NCR_SBRB/simplex/dng/affected_snvs.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_aff_dng = stats.itemfreq(aff_dng['snp'])\n",
    "unaff_dng = pd.read_table('/data/NCR_SBRB/simplex/dng/unaffected_snvs.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_unaff_dng = stats.itemfreq(unaff_dng['snp'])\n",
    "\n",
    "keep_me = [s for s, snp in enumerate(counts_aff_dng[:, 0]) if snp.find('chrX') < 0]\n",
    "counts_aff_dng = counts_aff_dng[keep_me, ]\n",
    "keep_me = [s for s, snp in enumerate(counts_unaff_dng[:, 0]) if snp.find('chrX') < 0]\n",
    "counts_unaff_dng = counts_unaff_dng[keep_me, ]\n",
    "\n",
    "my_max = np.max(counts_aff_dng[:, 1])\n",
    "max_str = counts_aff_dng[counts_aff_dng[:, 1] == my_max, 0][0]\n",
    "print 'Maximum frequency: %d in %d affected families' % (my_max, naff)\n",
    "print counts_aff_dng[counts_aff_dng[:, 1] == my_max]\n",
    "my_max = np.max(counts_unaff_dng[:, 1])\n",
    "print 'Maximum frequency: %d in %d unaffected families' % (my_max, nunaff)\n",
    "print counts_unaff_dng[counts_unaff_dng[:, 1] == my_max]\n",
    "print 'Counts %s shows up in unaffected: %d' % (max_str, counts_unaff_dng[counts_unaff_dng[:, 0] == max_str, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GATK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "cd /data/NCR_SBRB/simplex/gatk_refine\n",
    "\n",
    "# concatenate all affected trios and extract the SNPs\n",
    "rm affected_snvs.txt unaffected_snvs.txt;\n",
    "for f in `ls *_trio1_hiConfDeNovo.vcf`; do\n",
    "    cat $f | grep -v '#' - | awk 'BEGIN {FS=\"\\t\"; OFS=\":\"}; {print $1, $2}' - | sort | uniq >> affected_snvs.txt;\n",
    "done\n",
    "for f in `ls *_trio[2..4]_hiConfDeNovo.vcf`; do\n",
    "    cat $f | grep -v '#' - | awk 'BEGIN {FS=\"\\t\"; OFS=\":\"}; {print $1, $2}' - | sort | uniq >> unaffected_snvs.txt;\n",
    "done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum frequency: 2 in 21 affected families\n",
      "[['chr17:72762902' 2]\n",
      " ['chr18:15271172' 2]\n",
      " ['chr1:40229504' 2]\n",
      " ['chr20:29637674' 2]\n",
      " ['chr20:29637691' 2]\n",
      " ['chr2:89105006' 2]\n",
      " ['chr6:29968761' 2]\n",
      " ['chr6:32548712' 2]\n",
      " ['chr6:32557647' 2]\n",
      " ['chr6:32725367' 2]\n",
      " ['chr7:154002420' 2]\n",
      " ['chr7:154467985' 2]]\n",
      "Maximum frequency: 2 in 29 unaffected families\n",
      "[['chr12:11214231' 2]\n",
      " ['chr12:11214232' 2]\n",
      " ['chr12:50745851' 2]\n",
      " ['chr12:50745893' 2]\n",
      " ['chr12:50745894' 2]\n",
      " ['chr12:99139159' 2]\n",
      " ['chr16:33940111' 2]\n",
      " ['chr16:33940122' 2]\n",
      " ['chr1:142813302' 2]\n",
      " ['chr1:148902738' 2]\n",
      " ['chr1:206566826' 2]\n",
      " ['chr20:29638202' 2]\n",
      " ['chr2:9546136' 2]\n",
      " ['chr6:32083111' 2]\n",
      " ['chrX:116025284' 2]]\n"
     ]
    }
   ],
   "source": [
    "naff = 21\n",
    "nunaff = 29\n",
    "aff_gatk = pd.read_table('/data/NCR_SBRB/simplex/gatk_refine/affected_snvs.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_aff_gatk = stats.itemfreq(aff_gatk['snp'])\n",
    "my_max = np.max(counts_aff_gatk[:, 1])\n",
    "print 'Maximum frequency: %d in %d affected families' % (my_max, naff)\n",
    "print counts_aff_gatk[counts_aff_gatk[:, 1] == my_max]\n",
    "\n",
    "unaff_gatk = pd.read_table('/data/NCR_SBRB/simplex/gatk_refine/unaffected_snvs.txt',\n",
    "                               header=None, names=['snp'])\n",
    "counts_unaff_gatk = stats.itemfreq(unaff_gatk['snp'])\n",
    "my_max = np.max(counts_unaff_gatk[:, 1])\n",
    "print 'Maximum frequency: %d in %d unaffected families' % (my_max, nunaff)\n",
    "print counts_unaff_gatk[counts_unaff_gatk[:, 1] == my_max]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Some numbers\n",
    "\n",
    "First, how many DNVs do we get in each trio, per tool?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TrioDenovo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10033_trio1_denovo_v2.vcf 1188\n",
      "10033_trio2_denovo_v2.vcf 2353\n",
      "10042_trio1_denovo_v2.vcf 2399\n",
      "10090_trio1_denovo_v2.vcf 1884\n",
      "10090_trio2_denovo_v2.vcf 1712\n",
      "10094_trio1_denovo_v2.vcf 1895\n",
      "10094_trio2_denovo_v2.vcf 1846\n",
      "10128_trio1_denovo_v2.vcf 2201\n",
      "10128_trio2_denovo_v2.vcf 2157\n",
      "10131_trio1_denovo_v2.vcf 2033\n",
      "10131_trio2_denovo_v2.vcf 1125\n",
      "10131_trio3_denovo_v2.vcf 1164\n",
      "10131_trio4_denovo_v2.vcf 1333\n",
      "10153_trio1_denovo_v2.vcf 1119\n",
      "10153_trio2_denovo_v2.vcf 1912\n",
      "10153_trio3_denovo_v2.vcf 1007\n",
      "10164_trio1_denovo_v2.vcf 2141\n",
      "10164_trio2_denovo_v2.vcf 1869\n",
      "10173_trio1_denovo_v2.vcf 1897\n",
      "10173_trio2_denovo_v2.vcf 1916\n",
      "10178_trio1_denovo_v2.vcf 2577\n",
      "10178_trio2_denovo_v2.vcf 2423\n",
      "10182_trio1_denovo_v2.vcf 2347\n",
      "10182_trio2_denovo_v2.vcf 1333\n",
      "10182_trio3_denovo_v2.vcf 1830\n",
      "10197_trio1_denovo_v2.vcf 1281\n",
      "10197_trio2_denovo_v2.vcf 1179\n",
      "10215_trio1_denovo_v2.vcf 2311\n",
      "10215_trio2_denovo_v2.vcf 2270\n",
      "10215_trio3_denovo_v2.vcf 1122\n",
      "10215_trio4_denovo_v2.vcf 1177\n",
      "10369_trio1_denovo_v2.vcf 600\n",
      "10369_trio2_denovo_v2.vcf 628\n",
      "10406_trio1_denovo_v2.vcf 1560\n",
      "10406_trio2_denovo_v2.vcf 853\n",
      "10406_trio3_denovo_v2.vcf 1673\n",
      "10448_trio1_denovo_v2.vcf 1870\n",
      "10448_trio2_denovo_v2.vcf 1737\n",
      "10459_trio2_denovo_v2.vcf 1914\n",
      "1892_trio1_denovo_v2.vcf 1016\n",
      "1892_trio2_denovo_v2.vcf 941\n",
      "1893_trio1_denovo_v2.vcf 2816\n",
      "1893_trio2_denovo_v2.vcf 2607\n",
      "1895_trio1_denovo_v2.vcf 2016\n",
      "1895_trio2_denovo_v2.vcf 2080\n",
      "1976_trio1_denovo_v2.vcf 2164\n",
      "1976_trio2_denovo_v2.vcf 2225\n",
      "1976_trio3_denovo_v2.vcf 1126\n",
      "855_trio1_denovo_v2.vcf 1324\n",
      "855_trio2_denovo_v2.vcf 1940\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "cd /data/NCR_SBRB/simplex/triodenovo\n",
    "for f in `ls *_trio?_denovo_v2.vcf`; do\n",
    "    echo $f `grep -v '#' ${f} | sort | uniq | wc -l`;\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DenovoGear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10033_trio1_dnm.vcf 1417\n",
      "10033_trio2_dnm.vcf 2305\n",
      "10042_trio1_dnm.vcf 2234\n",
      "10090_trio1_dnm.vcf 1896\n",
      "10090_trio2_dnm.vcf 1810\n",
      "10094_trio1_dnm.vcf 1962\n",
      "10094_trio2_dnm.vcf 1845\n",
      "10128_trio1_dnm.vcf 1574\n",
      "10128_trio2_dnm.vcf 1896\n",
      "10131_trio1_dnm.vcf 2078\n",
      "10131_trio2_dnm.vcf 1538\n",
      "10131_trio3_dnm.vcf 1398\n",
      "10131_trio4_dnm.vcf 1663\n",
      "10153_trio1_dnm.vcf 1399\n",
      "10153_trio2_dnm.vcf 1710\n",
      "10153_trio3_dnm.vcf 1284\n",
      "10164_trio1_dnm.vcf 1909\n",
      "10164_trio2_dnm.vcf 1724\n",
      "10173_trio1_dnm.vcf 1825\n",
      "10173_trio2_dnm.vcf 1836\n",
      "10178_trio1_dnm.vcf 2310\n",
      "10178_trio2_dnm.vcf 2343\n",
      "10182_trio1_dnm.vcf 2157\n",
      "10182_trio2_dnm.vcf 1620\n",
      "10182_trio3_dnm.vcf 1775\n",
      "10197_trio1_dnm.vcf 1634\n",
      "10197_trio2_dnm.vcf 1414\n",
      "10215_trio1_dnm.vcf 1912\n",
      "10215_trio2_dnm.vcf 1900\n",
      "10215_trio3_dnm.vcf 1311\n",
      "10215_trio4_dnm.vcf 1262\n",
      "10369_trio1_dnm.vcf 636\n",
      "10369_trio2_dnm.vcf 635\n",
      "10406_trio1_dnm.vcf 1674\n",
      "10406_trio2_dnm.vcf 1187\n",
      "10406_trio3_dnm.vcf 1870\n",
      "10448_trio1_dnm.vcf 1824\n",
      "10448_trio2_dnm.vcf 1681\n",
      "10459_trio2_dnm.vcf 1824\n",
      "1892_trio1_dnm.vcf 1347\n",
      "1892_trio2_dnm.vcf 1213\n",
      "1893_trio1_dnm.vcf 2415\n",
      "1893_trio2_dnm.vcf 2475\n",
      "1895_trio1_dnm.vcf 2205\n",
      "1895_trio2_dnm.vcf 2255\n",
      "1976_trio1_dnm.vcf 2133\n",
      "1976_trio2_dnm.vcf 2328\n",
      "1976_trio3_dnm.vcf 1391\n",
      "855_trio1_dnm.vcf 1677\n",
      "855_trio2_dnm.vcf 1909\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "cd /data/NCR_SBRB/simplex/dng\n",
    "for f in `ls *_trio?_dnm.vcf`; do\n",
    "    echo $f `grep -v '#' ${f} | sort | uniq | wc -l`;\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GATK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10033_trio1_hiConfDeNovo.vcf 69\n",
      "10033_trio2_hiConfDeNovo.vcf 99\n",
      "10042_trio1_hiConfDeNovo.vcf 81\n",
      "10090_trio1_hiConfDeNovo.vcf 89\n",
      "10090_trio2_hiConfDeNovo.vcf 63\n",
      "10094_trio1_hiConfDeNovo.vcf 80\n",
      "10094_trio2_hiConfDeNovo.vcf 55\n",
      "10128_trio1_hiConfDeNovo.vcf 54\n",
      "10128_trio2_hiConfDeNovo.vcf 52\n",
      "10131_trio1_hiConfDeNovo.vcf 48\n",
      "10131_trio2_hiConfDeNovo.vcf 65\n",
      "10131_trio3_hiConfDeNovo.vcf 39\n",
      "10131_trio4_hiConfDeNovo.vcf 86\n",
      "10153_trio1_hiConfDeNovo.vcf 70\n",
      "10153_trio2_hiConfDeNovo.vcf 75\n",
      "10153_trio3_hiConfDeNovo.vcf 69\n",
      "10164_trio1_hiConfDeNovo.vcf 122\n",
      "10164_trio2_hiConfDeNovo.vcf 103\n",
      "10173_trio1_hiConfDeNovo.vcf 101\n",
      "10173_trio2_hiConfDeNovo.vcf 112\n",
      "10178_trio1_hiConfDeNovo.vcf 519\n",
      "10178_trio2_hiConfDeNovo.vcf 354\n",
      "10182_trio1_hiConfDeNovo.vcf 90\n",
      "10182_trio2_hiConfDeNovo.vcf 77\n",
      "10182_trio3_hiConfDeNovo.vcf 72\n",
      "10197_trio1_hiConfDeNovo.vcf 416\n",
      "10197_trio2_hiConfDeNovo.vcf 244\n",
      "10215_trio1_hiConfDeNovo.vcf 201\n",
      "10215_trio2_hiConfDeNovo.vcf 128\n",
      "10215_trio3_hiConfDeNovo.vcf 2023\n",
      "10215_trio4_hiConfDeNovo.vcf 186\n",
      "10369_trio1_hiConfDeNovo.vcf 948\n",
      "10369_trio2_hiConfDeNovo.vcf 996\n",
      "10406_trio1_hiConfDeNovo.vcf 119\n",
      "10406_trio2_hiConfDeNovo.vcf 81\n",
      "10406_trio3_hiConfDeNovo.vcf 103\n",
      "10448_trio1_hiConfDeNovo.vcf 69\n",
      "10448_trio2_hiConfDeNovo.vcf 70\n",
      "10459_trio2_hiConfDeNovo.vcf 93\n",
      "1892_trio1_hiConfDeNovo.vcf 91\n",
      "1892_trio2_hiConfDeNovo.vcf 61\n",
      "1893_trio1_hiConfDeNovo.vcf 95\n",
      "1893_trio2_hiConfDeNovo.vcf 105\n",
      "1895_trio1_hiConfDeNovo.vcf 82\n",
      "1895_trio2_hiConfDeNovo.vcf 74\n",
      "1976_trio1_hiConfDeNovo.vcf 91\n",
      "1976_trio2_hiConfDeNovo.vcf 96\n",
      "1976_trio3_hiConfDeNovo.vcf 69\n",
      "855_trio1_hiConfDeNovo.vcf 79\n",
      "855_trio2_hiConfDeNovo.vcf 104\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "cd /data/NCR_SBRB/simplex/gatk_refine\n",
    "for f in `ls *_trio?_hiConfDeNovo.vcf`; do\n",
    "    echo $f `grep -v '#' ${f} | sort | uniq | wc -l`;\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many of those are only show up in the affected trio?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TrioDenovo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10033 956\n",
      "10090 1238\n",
      "10094 1052\n",
      "10128 1080\n",
      "10131 1707\n",
      "10153 760\n",
      "10164 1515\n",
      "10173 1105\n",
      "10178 1720\n",
      "10182 1105\n",
      "10197 1078\n",
      "10215 960\n",
      "10369 402\n",
      "10406 826\n",
      "10448 1164\n",
      "1892 823\n",
      "1893 1625\n",
      "1895 1209\n",
      "1976 1150\n",
      "855 1075\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "suffix=denovo_v2.vcf \n",
    "cd /data/NCR_SBRB/simplex/triodenovo\n",
    "rm interesting_snvs_${suffix}.txt\n",
    "# figure out all family IDs\n",
    "ls -1 *_trio1_${suffix} > famids.txt\n",
    "sed -i -e \"s/_trio1_${suffix}//g\" famids.txt\n",
    "\n",
    "# for each family ID\n",
    "while read fam; do\n",
    "  # figure out how many trios we have\n",
    "   ntrios=`ls -1 ${fam}_trio?_${suffix} | wc -l`;\n",
    "   ntrios=$(($ntrios))\n",
    "   # if we have more than one (assuming the first one is affected)\n",
    "   if [ $ntrios -gt 1 ]; then\n",
    "      # get all SNVs in the affected trio in the family\n",
    "      cut -f 1,2 ${fam}_trio1_${suffix} | grep -v '#' - | sort | uniq > ${fam}_possible_snvs_${suffix}.txt;\n",
    "      # combine the vcf files of all unnafected trios\n",
    "      cat ${fam}_trio[2..$ntrios]_${suffix} > ${fam}_control_snvs_${suffix}.txt;\n",
    "      # for each possible SNV in affected trio, mark it as interesting if it's not\n",
    "      # in the unnafected trios\n",
    "      while read snv; do\n",
    "         if ! grep -q \"$snv\" ${fam}_control_snvs_${suffix}.txt; then\n",
    "            echo $snv >> interesting_snvs_${suffix}.txt;\n",
    "         fi;\n",
    "      done < ${fam}_possible_snvs_${suffix}.txt;\n",
    "      echo $fam `cat interesting_snvs_${suffix}.txt | wc -l`;\n",
    "      rm interesting_snvs_${suffix}.txt\n",
    "   fi;\n",
    "done < famids.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DenovoGear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10033 1128\n",
      "10090 1403\n",
      "10094 1349\n",
      "10128 1171\n",
      "10131 1614\n",
      "10153 922\n",
      "10164 1539\n",
      "10173 1226\n",
      "10178 1866\n",
      "10182 1165\n",
      "10197 1349\n",
      "10215 1042\n",
      "10369 431\n",
      "10406 947\n",
      "10448 1355\n",
      "1892 1074\n",
      "1893 1743\n",
      "1895 1499\n",
      "1976 1298\n",
      "855 1348\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "suffix=dnm.vcf \n",
    "cd /data/NCR_SBRB/simplex/dng\n",
    "rm interesting_snvs_${suffix}.txt\n",
    "# figure out all family IDs\n",
    "ls -1 *_trio1_${suffix} > famids.txt\n",
    "sed -i -e \"s/_trio1_${suffix}//g\" famids.txt\n",
    "\n",
    "# for each family ID\n",
    "while read fam; do\n",
    "  # figure out how many trios we have\n",
    "   ntrios=`ls -1 ${fam}_trio?_${suffix} | wc -l`;\n",
    "   ntrios=$(($ntrios))\n",
    "   # if we have more than one (assuming the first one is affected)\n",
    "   if [ $ntrios -gt 1 ]; then\n",
    "      # get all SNVs in the affected trio in the family\n",
    "      cut -f 1,2 ${fam}_trio1_${suffix} | grep -v '#' - | sort | uniq > ${fam}_possible_snvs_${suffix}.txt;\n",
    "      # combine the vcf files of all unnafected trios\n",
    "      cat ${fam}_trio[2..$ntrios]_${suffix} > ${fam}_control_snvs_${suffix}.txt;\n",
    "      # for each possible SNV in affected trio, mark it as interesting if it's not\n",
    "      # in the unnafected trios\n",
    "      while read snv; do\n",
    "         if ! grep -q \"$snv\" ${fam}_control_snvs_${suffix}.txt; then\n",
    "            echo $snv >> interesting_snvs_${suffix}.txt;\n",
    "         fi;\n",
    "      done < ${fam}_possible_snvs_${suffix}.txt;\n",
    "      echo $fam `cat interesting_snvs_${suffix}.txt | wc -l`;\n",
    "      rm interesting_snvs_${suffix}.txt\n",
    "   fi;\n",
    "done < famids.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GATK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10033 68\n",
      "10090 87\n",
      "10094 79\n",
      "10128 52\n",
      "10131 45\n",
      "10153 66\n",
      "10164 119\n",
      "10173 101\n",
      "10178 516\n",
      "10182 85\n",
      "10197 415\n",
      "10215 200\n",
      "10369 946\n",
      "10406 118\n",
      "10448 69\n",
      "1892 91\n",
      "1893 95\n",
      "1895 80\n",
      "1976 89\n",
      "855 77\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "suffix=hiConfDeNovo.vcf \n",
    "cd /data/NCR_SBRB/simplex/gatk_refine\n",
    "rm interesting_snvs_${suffix}.txt\n",
    "# figure out all family IDs\n",
    "ls -1 *_trio1_${suffix} > famids.txt\n",
    "sed -i -e \"s/_trio1_${suffix}//g\" famids.txt\n",
    "\n",
    "# for each family ID\n",
    "while read fam; do\n",
    "  # figure out how many trios we have\n",
    "   ntrios=`ls -1 ${fam}_trio?_${suffix} | wc -l`;\n",
    "   ntrios=$(($ntrios))\n",
    "   # if we have more than one (assuming the first one is affected)\n",
    "   if [ $ntrios -gt 1 ]; then\n",
    "      # get all SNVs in the affected trio in the family\n",
    "      cut -f 1,2 ${fam}_trio1_${suffix} | grep -v '#' - | sort | uniq > ${fam}_possible_snvs_${suffix}.txt;\n",
    "      # combine the vcf files of all unnafected trios\n",
    "      cat ${fam}_trio[2..$ntrios]_${suffix} > ${fam}_control_snvs_${suffix}.txt;\n",
    "      # for each possible SNV in affected trio, mark it as interesting if it's not\n",
    "      # in the unnafected trios\n",
    "      while read snv; do\n",
    "         if ! grep -q \"$snv\" ${fam}_control_snvs_${suffix}.txt; then\n",
    "            echo $snv >> interesting_snvs_${suffix}.txt;\n",
    "         fi;\n",
    "      done < ${fam}_possible_snvs_${suffix}.txt;\n",
    "      echo $fam `cat interesting_snvs_${suffix}.txt | wc -l`;\n",
    "      rm interesting_snvs_${suffix}.txt\n",
    "   fi;\n",
    "done < famids.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now find a few intersections. First, for all called SNPs:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, intersections for the SNPs in affected sibs only:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* RUN GATK WITHOUT FILTER? USE BOTH HICONF AND LOCONF DENOVOS!\n",
    "* GET NUMBER OF POSSIBLE DE NOVO PER FAMILY\n",
    "* MAKE VENN DIAGRAMS?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py2.7.10]",
   "language": "python",
   "name": "conda-env-py2.7.10-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
